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Abstract 

Propotional paths as summed up by the Goldscheider Rule (gr) - stating that given a constant 
strain rate, the evolution of the stress maintains the ratios of its components - is a characteristics 
of elasto-plastic motion in granular media. Barodesy, a constitutive relation proposed recently by 
Kolymbas, is a model that, with GR as input, successfully accounts for data from soil mechani- 
cal experiments. Granular solid hydrodynamics (gsh), a theory derived from general principles of 
physics and two assumptions about the basic behavior of granular media, is constructed to quali- 
tatively account for a wide range of observation - from elastic waves over elasto-plastic motion to 
rapid dense flow. In this paper, showing the close resemblance of results from Barodesy and GSH, 
we further validate GSH and provide an understanding for GR. 

PACS numbers: 81.40.Lm, 46.05. +b, 83.60.La 
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I. INTRODUCTION 

One focus of soil mechanical experiments is the stress evolution aij (t) for given strain rate 
Vij = |(Vjfj + Vjfi) and density p Three striking characteristics being observed at slow, 
elasto-plastic rates are (1) rate-independence, (2) the existence of a critical state i2|, and 
(3) proportional paths as summed up by the Goldscheider Rule (gr) [3[ . Rate-independence 
means that if the given strain rate Vij is a constant, the stress Cijlt) is a function of the 
strain Ey = f Vijdt = Vijt, and does not depend on the actual rate. 

The critical state is an expression of "ideal plasticity." Starting from an isotropic stress 
aij = PSij, and applying a constant shear rate v*, (* denotes the traceless part) - while 
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maintaining a vanishing trace vee = to keep the density p constant - a granular system 
will always go into an asymptotic, stationary state, in which the stress a'^j no longer changes 
with time, although Vij goes on providing a constant rate of deformation. We shall call 
this asymptotic state - characterized by the direction of the rate v*Jvs (where = v*jV*j) 
and the density p - the critical state. (The asymptotic state is more typically arrived at 
for given shear rate v*j, at constant pressure P = ^au or one of the principle stress value 
(Til, rather than the density. And there are some in the engineering community who insist 
on restricting the term critical state to the results of this second type of approaches. The 
narrower definition would be sensible if the respective asymptotic states were different. We 
do not believe this to be the case, for rather basic reasons, as will become clear in section HTl) 
The Goldscheider Rule or GR is a generalization of the critical state. First, it states that 
a granular system will converge onto the same critical state associated with v*j/vs and p, 
starting from any initial stress, not only an isotropic one. And second, it postulates the 
existence of asymptotic states also for cases of changing p - a point that we believe may be 
understood as follows: In the principal strain axes (^ii, £22, £^33), a constant v*j means the 
system moves with a constant rate along its direction, £11/6:22 = '^ii/'i^22) ^22/^^33 = ^22/^^3- 
This circumstance is referred to as a proportional strain path. In the stress space, a'^j is 
a stationary dot and does not move. Now, adding a constant Vu to the isochoric strain 
path, [v*j + t, we need to keep vu small compared to Vs, such that an initial stress 

has sufficient time to converge without breaching the random closest or loosest packing - 
the grains get crushed in the first case, and loose contacts with one another in the second. 
Then one expects the asymptotic state to be approximately given by the critical stress c^jlp) 
associated with the same isochoric strain path v*j t. As time passes, the density p will change, 
so will the critical state (jfj{p). But it will remain associated with v*jt. Interestingly, GR 
states that this stress path is also proportional, meaning crii{t) / (J22{t) , cr22{t) / o'ssit) also 
remain constant. 

The statements of GR as rendered by Kolymbas {4] are: (1) Proportional strain paths 
starting from the stress aij = are associated with proportional stress paths. (2) Propor- 
tional strain paths starting from aij 7^ lead asymptotically to the corresponding propor- 
tional stress paths obtained when starting at aij = 0. (A caveat: Although Goldscheider 
is a prudent and reliable experimenter, his data base is rather small (3|.) The initial value 
aij = is a mathematical idealization, neither easily realized nor part of the empirical data 
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that went into GR. So we take its core statement as: Applying a proportional strain path, 
there is an asymptotic line to which all other stresses converge. This stress path is also 
proportional, ^uch that proportional strain and stress paths form pairs. 

Barodesy [4] is a recent, impressively realistic constitutive model by Kolymbas, the orig- 
inator of hypoplasticity His purpose is to further improve the quantitative account 
of granular media's motion, and to reduce the considerable liberty when constructing hy- 
poplastic models. He did this by starting explicitly from GR. In the present paper, we 
take Barodesy (as specified below, in Sec IIII|) . as a constitutive equation with a reduced set 
of variables, which reflects highly condensed and intelligently organized empirical data, to 
which the results of GSH are compared. 

GSH is a theory of continuum mechanics derived from two notions that we hold to be 
the basic physics of granular media jo]. When constructing the theory employing general 
principles, we had little experimental data in mind, and certainly never needed to choose a 
subset of these. Therefore, if not totally wrong, GSH should be adequate and correct in a 
broad-ranged fashion. Until now, GSH has shown itself capable of accounting for phenomena 

nfi r n 

as diverse as static stress distribution [THSj , incremental stress-strain relatio n llOj , yield [UJ, 
12| . propagation and damping of elastic waves [l3|, elasto-plastic motion [ij], the critical 



state 



15| . shear band and fast dense flow 16|. Comparison to Barodesy is a further hard 



test for GSH, especially because any agreement could not possibly have been planned for. 
Moreover, GSH provides an understanding for GR, and embeds it among the many granular 
phenomena already understood within the framework of GSH. 

Finally, we amplify on the point why comparing GSH successfully to Barodesy validates 
the former. First, we note the qualitative difference between a physicists theory and an 
engineering one, which are in fact constructs of different raison d'etre: Physicists aim to 
first of all gain a qualitative understanding of a given system, while engineers want primarily 
to organize and mathematically condense experimental data gained in that system. When 
constructing a theory, physicists typically start from some notions about the basic behavior 
of a system, calling them motivation. If broadly validated, physicists will conclude these 
notions are right, considering this a gain in understanding. Two theories with different 
notions will contradict each other eventually, even if they initially agree with respect to 
some experiments. Two engineering theories may look very different mathematically, but 
if essentially the same set of data was used in constructing these theories, they will not 



4 



contradict each other starkly - though small deviations will generally exist. An agreement 
between GSH and Barodesy, showing that a theory deduced from some notions, hence quite 
possibly wildly off the mark, produces results that are seen in experiments, is therefore 
indeed a validation. On the other hand, there are always shades of gray - engineers with 
a fundamentalist's heart and physicists with a most pragmatic mindset. But all will agree 
that a mature theory should be both realistic and derived from the specific physics of the 
system under focus. 

We shall in future compare GSH to more constitutive models and experiments on elasto- 
plastic motion. Recent experiments include uniaxial tests [17], the settlement 18|, and 



systematic triaxial measurements [19j, though we do not expect either models or experi- 
ments to deviate strongly from Barodesy or each other, for the reasons mentioned above. 
A recent paper j2o| stands out, because it connects the constitutive model to particle- level 
properties. Generally speaking, one needs to heed the cautionary words by Schwedes |21| . 
that these setups are typically designed for engineering purpose, and the data may depend 
on operational details (such as skill of the experimenters and their level of training). So care 
has to be taken when adopting their data. A major difficulty, we believe, comes from the 
influence of the initial state, from the lack of information on boundary conditions, and from 
the presence of water (that GSH not yet considers). 

DEM simulations are nowadays a popular approach employed by both physicists and 



engineers, see eg. 22|-|32|. The best are often qualitatively perfect, but numerically different 
from experimental results. Therefore, a comparison to GSH requires us to treat them as 
different systems, with their own values for energy and transport coefficients. Unfortunately, 
the associated calibration process is time-consuming and laborious. 



II. THE EQUATIONS OF GSH 
A. Two-Stage Irreversibility 

The essence of granular physics, we contend, is encapsulated by two notions: two-stage 
irreversibility and variable transient elasticity. The first is related to the three spatial scales 
of any granular media: (a) the macroscopic, (b) the intergranular, and (c) the inner granular. 
Dividing all degrees of freedom into these three categories, we treat those of (a) differently 



5 



from (b,c). Macroscopic degrees of freedom, such as the slowly varying stress or velocity 
fields, are specified and employed as explicit state variables, but intergranular and inner 
granular degrees are treated summarily: Instead of being specified, only their contribution 
to the energy is considered and taken, respectively, as granular and true heat. So we do not 
account for the motion of a jiggling grain, only include its strongly fluctuating kinetic and 
elastic energy as contributions to the granular heat / Tgdsg, characterized by the granular 
entropy Sg and temperature Tg. Similarly, a phonon, or any elastic vibration within the grain, 
are taken as part of true heat, J Tds. Clearly, there are only a handful of macroscopic degrees 
of freedom (a), innumerable intergranular ones (b), and yet many orders of magnitude more 
inner granular ones (c). So the statistical tendency to equally distribute the energy among 
all degrees of freedom implies that the energy decays from (a) to (b,c), and from (b) to (c), 
never backwards. This is what we call two-stage irreversibility. 

Accounting for higher densities, when enduring contacts abound and granular jiggling is 
small, we expand wt to obtain wt = with Tg = dwr/dsg = Sg/b and wt ~ T^. There 

is no linear term because Sg,Tg = is an energy minimum. (The usual granular temperature 
Tg, defined as 2/3 of a grain's average kinetic energy, is useful only in the dilute limit, when 
the fluctuating elastic energy may be neglected, and wt ^ Tq ^ Tg.) Neglecting nonuniform 
situations, the balance equation for Sg reads 

TgdtSg = T]gVl + igVj^ " , (l) 

with vl = v*jV*j, and v*j the traceless part of Vij. Also: dt = The first two term on the 
right side accounts for viscous heating, the third for the leak of granular heat from (b) to 
(c). The viscosities ?7g,^g and the relaxation rate 7 are parameters of GSH and functions 
especially of Tg, p. We take joj 

Tg = Sg/ (pb), rig = riiTg, 7 = 7o + Ti^g, (2) 

noting that for what we call the hypoplastic regime of slightly elevated Tg, in which hy- 
poplasticity and Barodesy hold, 70 <^ liTg may be neglected. And we neglect ^gvj^, because 
rjgv"^ ^ ^gvjf^ in all typical experiments. For constant shear rate v*j, Eq ([1]) is a relaxation 
equation, with Tg quickly settling into its stationary value, 

Tg = ^/vi/li Vs- (3) 
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It is then no longer independent. The coefficients 1/6,771,71 are functions of p. 



b ~ (p, 



P) 



0.1 



71, ^1 ~ {Pcp- p) 



-1 



(4) 



cp 



where pcp is the closest packed density. We stand behind the Tg— dependence with much 
more confidence than that of the density, for two reasons: First, there are no comparably 
general arguments to extract the p dependence. Second, probably because of this, the 
observed dependence is not universal. The above dependence fits glass beads data, while 
7i ~ (Pcp — p) Vi ~ (Pcp ~ p)^^'^ seem more suitable for polystyrene beads, see jig . 

B. Variable Transient Elasticity 

Our second notion, variable transient elasticity, addresses granular elasticity and plastic- 
ity. The free surface of a granular system at rest can be inclined, or tilted. When perturbed, 
when the grains jiggle and Tg 7^ 0, the inclination will be reduced until the surface is hori- 
zontal. The stronger the grains jiggle, the faster this process is. We take this as indicative 
of a system that is elastic for Tg = 0, turning transiently elastic for Tg 7^ 0, with a stress 
relaxation rate that grows with Tg. A relaxing stress is typical of any viscous-elastic system 
such as polymers. The unique circumstance here is that the relaxation rate is not a material 
constant, but a function of the state variable Tg. As we shall see, it is this dynamically 
controlled, variable transient elasticity - a simple fact at heart - that underlies the complex 
behavior of granular plasticity. Realizing it yields a most economic way to capture granular 
rheology. 

Employing a strain field rather than the stress as a state variable usually yields a simpler 
description, because the former is in essence a geometric quantity, while the latter contains 
material parameters such as stiffness. Yet one cannot use the standard strain field Sij as a 
granular state variable, because the relation between stress and Eij lacks uniqueness when 
the system is plastic. A number of engineering theories divide the strain into two fields, 
elastic and plastic, Eij = Uij + e^j\ with the first accounting for the reversible and second 
for the irreversible part. They then employ Eij and e-J'' as two independent strain fields to 
account for granular plasticity [ssl. 




We believe that one should, on the contrary, take the elastic strain Uij as the sole state 
variable, as there is a unique relation between Uij and the elastic stress aij - if both are 
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related via the elastic energy: Shearing a granular system, part of the strain goes into 
deforming the grains, changing their elastic energy. The rest is spent sliding and rolling the 
grains. Taking Uij as the portion that changes the energy and deforms the grains, the elastic 
energy w{uij) is by definition a function of Uij alone. And since an elastic stress crij{uij) 
only exists when the grains are deformed, it is also a function of Uij. Therefore, we employ 
Uij as the sole state variable, and discard both Sij and e[f . Doing so preserves many useful 
features of elasticity, especially the (so-called hyper-elastic) relation, 

atj = -dw{uij,p)/duij\p. (5) 

This is derived in 6| but easy to understand via an analogy. Driving up a snowy hill slowly, 
the car wheels will grip the ground part of the time, slipping otherwise. (We assume a slowly 
turning wheel and quickly changing, intermittent stick-slip behavior.) When the wheels do 
grip, the car moves upward and its gravitational energy w^^°'^ is increased. If we divide the 
wheel's rotation into a gripping (e) and a slipping (p) portion, 9 = 9^'^^ + 9^p\ we know we 
may ignore 9^p\ and compute the torque on the wheel as dw^'^"''" / d9^^\ same as Eq.(l5]). How 
much the wheel turns or slips, how large 9 or 6'^^^ are, is irrelevant for the torque. 

The functional dependence of the energy density w{uij) is an input in GSH, as it cannot 
be obtained from general principles. The one we propose, because we find it both simple 
and appropriate, is specified below, in Sec III D[ Once it is given, so is the elastic stress, for 
which there is therefore an explicit expression, in terms of the state variables Uij and p. 

The evolution equation for Uj-,-, as derived in {sl, may be divided into that of the trace 
A = —uu and the traceless part m*,-, 

d,u*^ = (l-«)t;*.-ATX„ (6) 
dtA = (1 - a)f« + Ciiu*ikV*ik - >^iTgA. (7) 

Comparing them to the fully elastic equation, dtUij = Vij, we realize a is a gear shift 
factor, as a higher rate is necessary to achieve the same deformation; while ai is a dilatancy 
factor, accounting for the granular phenomenon that a shear flow leads to compression or 
decompression. Both a and ai are off-diagonal Onsager coefficients that depend on Tg, 
though we may take them as constant in the present context. 

The two terms ~ Tg are relaxation terms, accounting for the loss of deformation A, u*j 
(and the associated loss of the stress) when Tg is finite. The relaxation rate grows with 
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Tg, and is typically about 3 times as large for u*j as A, hence A ~ 3Ai [ij]. This is how 
variable transient elasticity is mathematically encoded. Replacing Tg with the shear rate Vg 
for the stationary case of Eq.Q, we find the above two equations explicitly rate independent. 



Denoting XTg = Xyrii/^iiVg = Avg, ^iTg = AiVg, we take, as tentatively as before, 

A, Ai, a, ai ~ p^p - p, (8) 

assuming that the plastic phenomena of relaxation, softening and dilatancy are no longer 
operative at p = pcp- For a given shear rate Vij, Eqs ( pfTI) are relaxation equations. Denoting 
= u*jU*j, v1 = v*jV*j, the system will converge onto the rate independent asymptotic value 
(denoted by a superscript of 

A ' ~ Ai M^Ai Vs' 

and u*j\'^ = v*j{u1/vs) implying and Vij have the same orientation and principal axes. 
Note also for vu ^ f^, the density dependence of A'^/u^ ^ tti/Ai cancels. 

Given the elastic strain Mjjl'^ = — and the density p, the elastic stress is also 

given. For vu = 0, this is simply the ideally plastic, stationary, critical state. The elastic 
strain and the associated stress do not change with time, because the deformation rate and 
the relaxation cancel, dtUij ~ dtcJij = 0. Calculating an approach to this state, starting from 



an isotropic stress and keeping the ai constant, the resultant curves for q = a^—cxi = (Ts a/3/2 
and the void ratio e = Pg/p — 1 {pg'- grain's density), against the strain in triaxial tests 
(cylinder axis along 3), resemble a textbook illustration of the critical state, see Fig[Tl 

For Vfj/vg finite but small, neither the direction nor the magnitude of Uij and aij change 
much, as conjectured in the introduction, though the mean stress will grow or decrease with 
the density. To see how it changes, and why it follows a proportional path, we need the 
explicit expression for the stress, specified in the next section. Sec III D[ But we can already 
see the reason for GR's second statement: Starting from any initial value Uij ^ the 
deviation Uij—Uij\'^ will, according to Eqs. ([6]) and ([7]), vanish exponentially. Clearly, variable 
transient elasticity is what lies behind both the critical state and GR. 

Finally, we note that the Cauchy stress, or total stress, a*"* is softer than the elastic one, 

a*f = (1 - a)a.j, (10) 

with the same factor as in Eqs. ([6]) or ([7]). This is a result of the Onsager reciprocity rela- 
tion [34r], but it may less formally also be seen as the flip side of the gear shift factor: when 




FIG. 1: Seemingly a textbook illustration of the critical state, this is the result of a GSH calculation: 
Shear stress q and void ratio e against the strain £3 in triaxial tests (cylinder axis along 3) , at given 
0"! and strain rate es/t, for an initially dense (po = 0.97pcp) and loose (po = 0.95pcp) sample. 
Insets are critical stress and density p"^ against tJi; with = 7GPa, A = Q.QB, pip = 0.85/9cp, 
ai = 750(1 - p/pcp), a = 0.7, A = 2850(1 - p/pcp), A/Ai = 3.3, and in SI units: 71 = 6 • 10^ 
r/i = 0.15 (ie. A = 114, ai = 30 for p = 0.96pcp, as in pj]). 

a higher rate is needed to achieve a certain elastic deformation, energy conservation requires 
the restoring force to be smaller by the same factor. 

For very high rates - such as present in heap flow or avalanches, the seismic pressure 
Pt ^ Tg vl resulting from violent jiggling of the grains, and the viscous stress, rjgV*j = 
riiTgV*j, again of shear rate squared, become relevant and destroy the rate independence. 
We neglect both in this paper - though not the viscous granular heating of Eq.(|3]). 
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C. Yield versus the Critical State 



In a space spanned by stress components and density, there is a surface that divides two 
regions in any granular media, one in which the grains necessarily move, another in which 
they may be at rest. This surface may be referred to as the yield surface. Equivalently, we 
may take the yield surface as the divide between two regions, one in which elastic solutions 
are stable, and another in which they are not - clearly, the medium may be at rest for a 
given stress only if an appropriate elastic solution is stable. Since the elastic energy of any 
solution satisfying force equilibrium VjCTjj = is an extremum the energy is convex 
and minimal in the stable region, concave and maximal in the unstable one — in which 
infinitesimal perturbations suffice to destroy the state of rest. 

Yield is clearly a completely distinct concept from the critical state as discussed above - 
one a static phenomenon, the other dynamic. The first is a convexity transition of the elastic 
energy, to be probed by quasi-static motion at vanishing Tg, say by slowly tilting a plate 
with a layer of grains. The second is a stationary solution of the relaxation equation for 
the elastic strain Uij at given shear rate, and relevantly, at an elevated Tg. It is comparable 
to the stationary solution of any diffusion equation. The yield and critical shear stresses 
are frequently similar in magnitude, but the yield stress needs to be larger than the highest 
shear stress achieved during the approach to the critical state. Otherwise, the system will 
abandon the approach and develop shear bands instead, see Figj2] below. 

Many believe that the approach to the critical state is accomplished at low enough shear 
rates to be considered quasi-static. We contend that a quasi-static motion is one that visits 
a series of static, equilibrium states, with Tg — )■ 0. The rate of dissipation must be neghgibly 
small. The rate-independent, hypoplastic motion taking place during an approach to the 
critical state maintains an elevated Tg that allows irreversible, dissipative relaxation of the 
elastic strain. In the critical state, this dissipative process, having the same magnitude as 
the reactive (ie. elastic) deformation rate, certainly cannot be neglected. 

D. Granular Elastic Energy 

The elastic energy density w is a function of the three independent strain invariants, 
A, Us and Ut = \/u*^u\-u*^. For granular materials, the following expression is appropriate 
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in many respects, 

^ = SyA0A^ + l«^-||), (11) 

and was instrumental in achieving the agreements with all the granular phenomena men- 
tioned above, especially static stress distribution, incremental stress-strain relation, and 
elastic waves. Varying the coefficient x? the yield surface changes to resemble different yield 
laws, including Drucker-Prager, Lade-Duncan, Coulomb, and Matsuoka-Nakai, see [l^. For 
qualitative considerations, however, it is frequently sufficient to set x = 0. We then find the 
elastic energy convex only for 

M,/A<v^, or ajP<^i, (12) 



where P = o"fcfc/3, = \/of^Mf^- The energy w turns concave if this condition is violated. 

We keep a finite x fo^' the rest of this paper, taking it along with ^ ?^ 5/3 as density 
independent. But B{p) is specified as ^ 

B = Bo[ip-p)/ipc,-p)r^ (13) 

with p = (20p£p — llpcp)/9, and Bo > 0. This expression accomplishes three things: • The 
energy is concave for any density smaller than the random loose one p^p, implying no elastic 
solution exists there. • The energy is convex between the random loose density pip and the 
random close one pcp, ensuring the stability of any elastic solutions in this region. In addition, 
the density dependence of sound velocities as measured by Harding and Richart [35| is well 
rendered by \/B. • The elastic energy diverges, slowly, at pcp, approximating the observation 
that the system becomes orders of magnitude stiffer there. 
The elastic stress aij = —dw/duij may be written as 

and calculated employing Eq.l ITTi) . Using Eqs.l fT^ it can also be shown that for any isotropic 
energy the stress and elastic strain tensors have same principal directions. And since the 
critical elastic strain is colinear with the strain rate, all three have the same principal axes 
asymptotically. The stress eigenvalues o"i^2,3 are 



a,- = By/X 



1 x\u 



,2 



X <4 2 ( A\ ^ 3x ( . A^' 
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where Ui denote eigenvalues of Uij. From (fT5|) . the following relations between the triplet of 
strain invariants {A,Us,Ut) and stress invariants iP,as,o't = \/ '^ik'^lj'^ji ) holds: 




We are now in a position to understand that the proportional stress path is in fact a result of 
certain coefficients (or combinations of them) not depending on the density: The above three 
formulas show that P/cr^ and at/cTs depend on x, ^, A^/u'^, ul/u^. If all four are independent 
of the density, the stress path is proportional, with the stress magnitude growing with p, as 

stress magnitude ~ (1 - a)B{p) A^/^, (19) 

cf. Eqs fllOmsp . These four quantities are indeed density independent for vi,V2 ^ f«: First, 
because of ul/u^ = Vt/vs, see the first equation after Eq ([9]), or alternatively Eq ( l43l) below, 
the quantity ul/Ug depends only on the direction of the shear rate. Second, A'^/ul ^ 
tti/Ai have been taken above as density independent. 



III. THE BARODESY MODEL 



In soil mechanics, granular dynamics is frequently modeled employing the strategy of 
rational mechanics^ by postulating an algebraic function - of the stress cXjj, strain rate 
and density p - such that the constitutive relation, {dt + v^V kl^^ij = ^ijii^ij, Vke, p) holds. 
(Instead of the density, one can equivalently take the void ratio e = Pbuik/p — with p^uik 
the bulk density of the grains.) Together with the continuity equation dtp + Vipvi = 0, 
momentum conservation, dt{pvi) + Vj(crjj + pviVj) = 0, it forms a closed set of equations for 
p, aik and the velocity Vi. Barodesy as proposed by Kolymbas [4] is such a model. Note that, 
in comparison to GSH, the state variables Tg and Utj have been eliminated by taking certain 
limits, such as stationarity or rate-independence. It is therefore a constitutive equation with 
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a reduced number of variables. Barodesy is defined by the following expressions: 



^1 — C3 



(J, 



\C3/2 



mn ^ mn I 



CAVkk 



+ 



where the direction and magnitude of the asymptotic stress are given respectively as 

Vkk . . ( Vi 



(20) 



Tij — 



, --?>ij + ci exp C2 



1 + ec 
1 + eoc 



exp 



(l-C3)/2 



C4 (1 - C3) 



(21) 
(22) 



Similar to GSH, the asymptotic stress has the same principal axes as the strain rate 
ie. Tij is diagonal if % is. There are 6 dimensionless parameters, the values of which are: 



ci = -1.7637, C2 = -1.0249, C3 = 0.5517, 
C4 = -1.174, C5 = -3.26, eoc = 0.75, 



(23) 



in addition to one with dimension, cq, in Pa. 



IV. STRAIN PATHS OF CONSTANT DENSITY 



In this section, we evaluate GSH and Barodesy analytically. This is possible only for 
vu = and constant density. So we are dealing in fact with the critical state here, as 
discussed above a special case of GR. A proportional deformation path is given by all three 
eigenvalues of vij {t) remaining proportional. We decompose the strain rate tensor as 



^ vi ^ 

\0 V3 J 



(24) 



with Vi + V2 + V3 — 0, ■'^vf + v^ + v^ — 1, and v = y/vikVik > denoting the magnitude. 
Proportional paths imply the constancy of vi^2,z- Next, we employ the strain Lode angle, 

= (1/3) arcsin (y&vl/vl) , (25) 
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(where Vt = ^v*^vl-v*^, and Vs = v because Vkk = 0) to express f 1,2,3, 



V2 



V3 



sin — 



sinLp 



^sin(L. + £). 



(26) 
(27) 
(28) 



Similarly, we can write the stress tensor as 



^ ai ^ 



Cij = (tR 



(29) 



3^2 
\0 as J 

with R a rotation matrix - equal to the unit matrix for the asymptotic states in both GSH 
and barodesy. ^1,2,3 need to be expressed by two angles, because ai + 5^2 + ?3 7^ 0. We 
take them as the stress Lode angle L and the friction angle (, as defined in the Appendix, 
Eqs.(|A4]|A5|), 

cos ( 2 sin C . 



^3 



_ , _ sin ( L - ^ 

cos ( 2 sin ( 
—= — sm L, 

cosC 2sinC . , 7r\ 



(30) 
(31) 
(32) 



Proportional path implies time-independent L, (. [The relation between the angles L, ( and 
the stress invariants P,(Js-,(Jt are given in the Appendix, Eqs.f lA4|A5p .] The association 
between the strain and stress paths may be given as 

L = L{L,) C = C(^.) (33) 
We now calculate the stress evolution for proportional strain paths, 

^ sin (L^ - f ) ^ 

-sinL^ , (34) 

Y sin(L, + f)^ 

obtained by inserting ( I26|27f28p into (12^ . taking v = const. Both GSH and Barodesy deliver 
analytical expressions. 
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A. Results from GSH 



The GSH equations can be solved as follows. First, inserting the strain rate into 
Eq.dl]), and noting that v = Vg = const, we have that the solution for Tg is: 

Tg = v^rixhx tanh [vt ^Ti^i ) (35) 
= v\/r]x hi tanh {e a/7i?7i ) 

Here the initial condition Tg{t = 0) = is assumed. The notation e = j^vdi' = vt is 
magnitude of total strain. Clearly 1/ (vy/^jifji) (or 1/a/7i77i) is the time (or strain) scale 
needed for Tg going to its saturation (asymptotic) value. 

Inserting (l35l) into Eq.([6]), we obtain the deviatoric elastic strain, 

^ + (1 - ^s) V*, /o cosh"/^^ jvyimt') dt' 

cosh-^/^^ ^ 
where u*-q = u*- (t = 0) is the initial strain. Because coshx — )■ for x — oo, the initial 



strain mT q decays as e The dimensionless parameter 



A = X^/^iJTi (37) 

is typically ~ 114 for the intermediate void ratio of e = 0.65, see jo]. So the decay is fast, 
ending at about 1% strain magnitude e. Inserting ( 135|) into ([7]), we have 

A = Ao + J>y)coshY-C"V7;^'V''^ (33) 

cosh-^^/T^ (VW^) 
, »r;co^r;c + (!-«) v' /p cosh^/^^ (t; VTl^tQ 

ft = ai w . . ■ (39) 

cosh-"/^^ (VtI^^) 

So the initial bulk strain Aq decays with e~^^^, more slowly, as Ai A/3. For t oo, the 
strain ( l36|38l) becomes stationary: 



A ^ A. = ^lli^. (41) 
AAi ^ ^ 

The associated three invariants are 

< = ^' ^ = -> (42) 

(43) 
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1 — a 


_ Ai 


A ' 


A^ ai 


_ ^* _ 


sini/3 (3^^ 


f 


gl/6 



where the third expression is obtained with Eq. (l25|) . Inserting these into (fT5|) . we have the 
following principal stresses: 

3/2 



^1 



B (\ 
1 



a 



A 

+ 1^ sin (3L,) - 2Csm (l, - ^) + 3xC' sin (l, - ^ 



C 



1 



(44) 



^2 



H ^ sin (3Le) + 2C sin — 3yC^ sin 



+ ^ sin (3L,) - 2C sin (l, + ^) + 3xC^ sin (l^ + 
2v6 V 3/ V 



(45) 



(46) 



where C = Ai/ai. In the 7r-coordinates (defines as tti = (cr^/P) sin (L + |) , 7^2 = 
— ((Js/P) cos (-^^ + f ), see the Appendix for more details), the critical stress has more com- 
pact expressions: 



7r2 



6C2;^sin (2L, + f ) - 476^ cos (L, - f ) 

2 + 76^2 + sin 3L, 
6^2;^ cos (2L, + f) -4y6C sin (Lf. — ^) 



(47) 

^ ^ • (48) 

2^6^ + 76^2 + ;^C3 sin 3L, 

When the Lode angle varies from to 27r, the loci given by Eqs.( H7|48|) give a triangle-like 
curve, as shown by the full line in Figj2l The curve is determined by the three parameters: 
C, ^, X; and reduces to a circle if x = 0. 



B. Results from Barodesy 

The critical surface of the Barodesy is obtained by taking dtaij = and Vkk = 0. In this 
case, Eq.(l20l) reduces to e = Cc and 



(Tij 



exp ( . (49) 



exp / ■ exp ■ , 
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FIG. 2: Loci of asymptotic states in the vr-plane, as obtained, respectively, from GSH (full) and 
barodesy (dots). Dashed curve is the static yield surface, ie. the convexity surface of the energy 
w, Eq.pT]). 



Inserting the strain rate Eg. (1341) in the principal axes, we have 



0-2 



0^3 



(0-10-20-3)^^^ exp 



0-1 0-20-3)^^^ exp I -C2W^sinL 



02^1 sin (l,- 



((Ti(T20-3)^/%Xp 



C2\/^sin (^Le + !j 



or in the vr-coordinates: 

3 



vr2 



exp [a/2 |c2| cos (Le - 1)] 



\/2 exp (-\/2 IC2I cos Lg) + exp [a/2 |c2| cos (L^ — f )] +1 



3 2 exp (a/2 


|C2 


COS Le) — exp [v^ 


|C2 


|C0S (Le-f)] -1 


\/6 exp (\/2 


C2I 


cos Le) + exp [a/2 


C2I 


COS (Le-f)] +1 



(50) 
(51) 
(52) 



In barodesy, the critical surface is determined by the parameter C2, and is triangle- 
the vr-plane, see Figj2l which is the same curve as in Fig.l of the second reference of 



(53) 

(54) 

ike in 
4|. 



Transforming both the GSH expressions of Eqs. (l47ll48l) and the Barodesy ones of 
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FIG. 3: Angle association of Eq.(33), calculated employing GSH (full) Barodesy (dots). 

Eq.( l53|54|) into the angles L^C, using the formula given in the appendix, we retrieve the 
association of Eq.f l33p . as shown in Figj3l again with great similarity between both theories. 

In contrast to the last two figures that contain only asymptotic information, FigJH shows 
the evolution of three stress eigenvalues, starting from an initially isotropic stress state. 
The numeric calculation employs GSH, Eqs. (161I71IT]) and Barodesy, Eqs.(l20l). for = 15°. 
The transient behavior is clearly somewhat different, it contains an oscillation in GSH (full), 
but is monotonic in Barodesy (dashed). The discrepancy is probably due to the (correct) 
nonmonotonic behavior of the pressure and shear stress in GSH. 

Although the expressions from barodesy, Eqs.( l53]l54l) . and GSH, Eq. ( l47|48l) . are rather 
different, the relevant plots are not. Yet to achieve this agreement, hardly any fiddling with 
the parameters was necessary. The Barodesy parameters were simply taken from Eq.( l23|) ; 
the GSH parameters are essentially the same as we employed them before: ^, x are part of 
the energy and represent static parameters. We took ^ = 5/3 as we have mostly done before, 
and took x = 0.26. (In [l^, we took x = 0-2- This slight change perfected the agreement 



of Fig.l that we could not resist.) We also took C = AAi/ (Aai) = 1 



equivalently took A = Aa/^TiTti ~ 114, Ai ~ 3A/10, and ai ~ 33 



in [l4 



lere. Previously, we 

n 

115[, separately. 
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Stress evolutions for proportional strain path for A=15 
Initial stress: P =0.4 kPa a =0 
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FIG. 4: (a) Evolution of the stress eigenvalues (normalized by the pressure P) with the total strain 
e, as obtained with GSH (full) and Barodesy (dashed), (b) Evolution of P (normalized by the initial 
pressure Pq) with e. The Barodesy curve is monotonic, The GSH one is not. This explains the 
difference in the transient regime. 

V. STRAIN PATHS WITH DENSITY CHANGE 



If the strain path contains a small vu-, the density and void ratio will change, as will the 
magnitude of the stress, according to Eq. f|T9|) in GSH, and to Eq. fl22|) in barodesy. Again, in 
spite of the different expressions, the curves are similar, at least qualitatively, see Fig.([S]). 
The convergence onto the asymptotic state is depicted in Fig.(l6l). Following Kolymbas' 
papers, we have also computed 4 figures each for (drained) triaxial and oedometric tests 
employing GSH, see Figl7]and[8l 
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FIG. 5: The relation between the void ratio e and the stress magnitude a of the asymptotic state 
(normaUzed by cr at e = 0.65), as given respectively by GSH and barodesy. With vu 7^ 0, the two 
quantities e and a will change in tandem, according to either of these two curves. 

VI. CONCLUSION 

In comparing GSH to barodesy, we set out to achieve two goals: To validate GSH, and to 
provide a transparent, sound understanding for GR. Both goals were reached. GSH is vali- 
dated because it yields similar results for various key quantities as barodesy, achieving better 
agreement than would be reasonably expected, without much fiddling with parameters. 

Conversely, the understanding of Goldscheider Rules and Barodesy comes from the 
physics of GSH. The theory has (for the range of shear rates typical of soil mechanical 
experiments) three state variables and two constituent parts. The state variables are the 
density p, granular temperature Tg (quantifying granular jiggling), and the elastic strain 
tensor u (accounting for the coarse-grained deformation of the grains). The two constituent 
parts are first the explicit expression of the stress tensor o", a function of m, p that is obtained 
from the elastic energy; and second a rate-independent relaxation equation for it, derived 
from the notion of variable transient elasticity. Given any initial Uq, the system will always 
converge onto the stationary solution vf^ as prescribed by the relaxation equation, vf^ is a 
function of the constant strain rate, or equivalently, of the proportional strain path's di- 
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FIG. 6: Circumstances are the same as in Fig. 3, except that, first, v^k 7^ and the density varies 
linearly, as given in the inset of (b). Second, only the GSH curve is displayed, not the Barodesy one 
(because given Fig l4|5l they cannot be that different). Convergence of the stress ratios takes place 
quickly, from where on the stress path is proportional. The first part of the pressure change occurs 
during the convergence, the second part, where the pressure change associated with positive and 
negative Vkk diverge, belongs to the asymptotic state. 

rection, and may be identified as the asymptotic, critical state for isochoric paths, vu = 0. 
This convergence, a consequence of the relaxation equation, is closely related to variable 
transient elasticity, and hence a generic aspect of granular behavior. 
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0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 

FIG. 7: In the geometry of triaxial tests, various quantities are computed employing GSH, as 
functions of the strain £xx, holding axx = (^yy constant. (The axial direction is z. The case with 
an initially higher density is rendered in solid lines, the looser one in dashed lines.) These are: 
(a) deviatoric stress q = azz — ctxx', (b) void ratio e; (c) volumetric strain e^,; (d) the friction ange, 
siiKpm = q/ {2axx + q)- We chose: a,ai,X ~ (1 - p/pcp)^'^ and ??i,7i ~ (1 - p/pcp)~^. 

Given u'^ and the density, the stress is also fixed. Its form, however, depends on the 
expression for the elastic energy that is material dependent and less robust. If the strain 
path -0* t is isochoric, with vu = 0, the asymptotic stress state is a constant of time, but 
a function of u^, or equivalently, of the strain path's direction. As the path varies, the 
associated stress states lie within a triangle, as depicted in Figj2l 

If the shear rate is a sum of v* t and a small vu, the asymptotic state is (cum grano 
salis) still given by the u'^ associated with v* t, though the density will now change. The 
asymptotic stress is therefore a function of the same u'^ and a changing density, hence no 
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FIG. 8: Oedometric test starting with an initially loose density, again computed with GSH. The 
following four quantities are calculated as functions of the axial Gzz while holding the lateral strain 
^xx = £yy fixed: (a) of a axial stress-strain curve; (b) axx of a stress path; (c) e^z versus time; 
(d) void ratio e versus time. Note Ezz moving upwards implies a compaction of the system. 

longer a constant. That the stress path is also proportional, that only the magnitude of the 
stress changes with time, not the ratios of its eigenvalues, is the least robust part of GR, 
because it depends on the density dependence of certain coefficients canceling. 

Constructing a constitutive relation, specifying {dt + VkVk)c^ = ^{o-,v,p), is only for 
someone with vast experience with granular media and deep knowledge of how they behave. 
That GSH - derived from two simple notions of what the two basic elements of granular 
physics are - yields an equivalent account, is eye-opening, and the actually amazing fact of 
the presented agreement. 



24 



Acknowledgments 



Helpful discussion with, and critical reading of the manuscript by, Dimitris Kolymbas are 
gratefully acknowledged. 



Appendix A: Tensor decomposition 

A 3 X 3 symmetric tensor, e.g. the stress tensor be decomposed into two parts: 

a spatial rotation O and a part which is invariant under any rotation. In most analysis we 
are interested mainly in the three invariants. There are usually various ways to represent 
the invariant triplet. One of which is 

P = 'Jkk/^, (Al) 



'^s = V^k^ (A2) 



<^t ^ {l^lk^lfh- (A3) 



Another is {P,L,C,) where 

L = - arcsin ( \ , (A4) 
C = arctan ^— = arcsin ^ — arccos ^-s/S— ^ , (A5) 

are two angle variables {a = ■^JoikOik — a/ct^ + 3P^). In soil mechanics L is usually called 
the Lode angle of stress. The angle C can be interpreted as a "friction angle" (because 
it represents the ratio between shear force and pressure). Moreover we can also use the 
three eigenvalues (cxi, (T2, era) of the stress tensor as an invariant triplet, which are related to 

(P, at) by 

<T2 = P-H^sinL, (A7) 
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where L is given by (lA4p . In soil mechanics, it is also usual to define the two coordinates 
(7ri,7r2) in the so called vr-plane, 

Ti = (A9) 

2cri-cr2-0-3 / . , „x 

Inserting (lA6]IA7yA7p into ( 1A9]IAT0|) . we have 

TT^ = ^sin('L + ^) , (All) 



TTa = -^cos (l^^) . (A12) 



P V 6 
^cos (L+g 

With the help of Eqs. (]AltlA12|) we can readily transform among the invariant triplets: 
{P,as,(Tt), {P,L,Q, (c^i, cr2, CTs), TTl, 7r2). Similar decompositions apply for the elastic 
strain Uij, total strain , strain rate tensor Vij etc., only note that the first invariant is fre- 
quently defined with a factor different from that of P. 
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